
***********************************************************

* (1.0) Set up the analysis sample
***********************************************************

	use 	"$f_2PKbirthpanel", clear	

	keep if rural ==1	
	merge m:1 prov year using  "$d_data\IMR.dta", gen(IMRmerge)

	gen 	birth_year = .
	replace birth_year = year
		
	cap drop 	birth_month
	gen 		birth_month = .
	replace 	birth_month = pregmonth 
	gen 		marriage_year = .
	replace 	marriage_year = 1900 + year_marriage 
	gen 		marriage_age = age_marriage 

	cap drop order
	sort 	 	id birth_year birth_month
	by 			id: gen order = _n
	
	gen 		months_since_mar = bcmc-marcmc
	replace 	months_since_mar = birth_month + 12*(birth_year-(1900+year_marriage)) - month_marriage 
	
	keep if 	months_since_mar > 0 & months_since_mar != .
	
	cap drop 	temp
	gen 		temp = months_since_mar if order == 1
	bysort 		id: egen interval1 = mean(temp)
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval1 if order == 2
	bysort 		id: egen interval2 = mean(temp)
	replace 	interval2 = . if interval2 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval2 - interval1 if order == 3
	bysort 		id: egen interval3 = mean(temp)
	replace 	interval3 = . if interval3 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval3 - interval2 - interval1  if order == 4
	bysort 		id: egen interval4 = mean(temp)
	replace 	interval4 = . if interval4 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval4 - interval3 - interval2 - interval1  if order == 5
	bysort 		id: egen interval5 = mean(temp)
	replace 	interval5 = . if interval5 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval5 - interval4 - interval3 - interval2 - interval1  if order == 6
	bysort 		id: egen interval6 = mean(temp)
	replace 	interval6 = . if interval6 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval6 - interval5 - interval4 - interval3 - interval2 - interval1  if order == 7
	bysort 		id: egen interval7 = mean(temp)
	replace 	interval7 = . if interval7 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval7 - interval6 - interval5 - interval4 - interval3 - interval2 - interval1  if order == 8
	bysort 		id: egen interval8 = mean(temp)
	replace 	interval8 = . if interval8 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval8 - interval7 - interval6 - interval5 - interval4 - interval3 - interval2 - interval1  if order == 9
	bysort 		id: egen interval9 = mean(temp)
	replace 	interval9 = . if interval9 < 9
	
	cap drop 	temp
	gen 		temp = months_since_mar - interval9 - interval8 - interval7 - interval6 - interval5 - interval4 - interval3 - interval2 - interval1  if order == 10
	bysort 		id: egen interval10 = mean(temp)
	replace 	interval10 = . if interval10 < 9
	
	cap drop _merge
	merge 	m:1 prov year using `"$d_data\ProductionData.dta"'
	
	cap drop _merge
	merge m:1 prov using `"$d_data\sentdown_mean.dta"'
	drop _merge 
	
	gen 	culturalrev = 0 
	replace culturalrev = mean_sentdown if year >= 1966 & year <=1969
	
	do `"$d_do\LLF_StackData.do"'

	keep id order months_since_mar year_marriage month_marriage marcmc birth_year birth_month interval1 interval2 interval3 interval4 interval5 interval6 interval7 interval8 interval9 interval10 LLFyear prov m_year_birth year educ marriage_age U5MR_5yr_ave grain agric_production gdp culturalrev pop_total_agric FE_prov FE_t StackDataID

	drop if order>10
	
	gen eventyear = birth_year - LLFyear
	drop if eventyear>7
	drop if eventyear<-8
	char FE_prov[omit] 108 
	
	global controls1 	"i.educ  U5MR_5yr_ave"
	global controls2 	"grain agric_production gdp pop_total_agric"

***********************************************************
* (2.0) Logit model using quarterly survival time 
***********************************************************

	gen agefirstbirth = marriage_age + ((interval1*3)/12)			
	egen uniqueid = group(id StackDataID)
		
	replace year_marriage = 1900 + year_marriage

	gen years_since_LLF = year - LLFyear
	cap drop period
	gen period = 0 if years_since_LLF <= 0 
	replace period = 1 if years_since_LLF >= 1 & years_since_LLF <= 4
	replace period = 2 if years_since_LLF >= 5

	local i = 1 
	
	sum interval`i' 
		local overall = r(mean)
	sum interval`i' if period == 0 
		local pre = r(mean)
	sum interval`i' if period == 1 
		local early = r(mean)
	sum interval`i' if period == 2  
		local late = r(mean)

	mat sumstats = `overall', `pre', `early', `late' 
	
	forval i = 2/4 { 
	
	sum interval`i'  
		local overall = r(mean)
	sum interval`i' if period == 0 
		local pre = r(mean)
	sum interval`i' if period == 1 
		local early = r(mean)
	sum interval`i' if period == 2  
		local late = r(mean)
		
	mat temp = `overall', `pre', `early', `late' 
	mat sumstats = sumstats \ temp
	} 

	
	tempfile main
	save `main', replace 

	************************
	* (3.1) First births
	************************
	preserve
	
	keep if order == 1  & year_marriage >= 1960 & birth_year <= 1979
	*keep if agefirstbirth >= 24
	keep marriage_age m_year_birth educ U5MR_5yr_ave $controls2 birth_month birth_year LLFyear year uniqueid interval1 prov marcmc FE_prov FE_t 
	
	* I want to make this a woman-quarter level dataset so I can allow LLF to vary over the interval
	* Define quarterly survival (starting at 9months postpartum)
	local interval = "interval1" 
	gen 	interval0 = interval1 + (marriage_age*12) - (15*12)
	lab var interval0  "Interval in months from 15 to first birth" 
	
	local interval = "interval1" 

	* gen quarter indicators that capture survival time - they should be zeros until the quarter the interval ends, when it should be 1
	gen 	quarter1=0 if `interval'>=12 & interval1!=.
	replace quarter1=1 if `interval'<12
	
	forval i = 2/56 { 
		local j 	= 3*`i' 
		local low 	= 9 + `j' - 3
		local high 	= 9 + `j'
	
		gen 	quarter`i'=0 if `interval'>=`high' & interval1!=.
		replace quarter`i'=1 if `interval'>=`low'  & `interval' < `high'	
	}
	
	replace quarter56 = 1 if `interval' >= 177
	
	* Now can identify the year and month at the start of the interval 
	* Note that first month at risk is 9 months after the previous event, but the start is the start of the entire birth interval
	gen end_cmc = 12*(birth_year-1900) + birth_month
	gen start_cmc = end_cmc - interval1 
	
	forval i = 1/56 { 
		cap drop q_cmc`i'
		gen q_cmc`i' = .
		
		local q = `i'*3 - 3
		replace q_cmc`i' = start_cmc + 9 + `q' if quarter`i'  != .
		
		cap drop q_year`i'  
		gen q_year`i' = . 
		replace q_year`i' = q_cmc`i'/12 + 1900 if quarter`i'  != .
		
		cap drop q_LLF`i' 
		gen q_LLF`i' = . 
		replace q_LLF`i' = 0 if quarter`i' != . 
		replace q_LLF`i' = 1 if q_year`i' >= LLFyear & quarter`i' != . 
		} 
		
	drop q_cmc*  end_cmc start_cmc
	
	reshape long q_LLF q_year quarter, i(uniqueid) j(interval_q)
	replace q_year = floor(q_year)
	rename quarter birth
	rename interval_q quarter
	drop if birth == .
	
	stset quarter, failure(birth) id(uniqueid)
		
	tab quarter, gen(quarter_)
	forval i = 1/56 { 
		gen LLF_quarter_`i' = q_LLF*quarter_`i' 
		} 
		
	gen birth_year2LLF = birth_year-LLFyear
	
	keep if q_year >= 1964
	keep if birth_year2LLF >= -4	
	keep if birth_year2LLF <= 9		
	
	xi: logit 	birth q_LLF LLF_quarter_1-LLF_quarter_12 LLF_quarter_14-LLF_quarter_56 quarter_1-quarter_12 quarter_14-quarter_56 ///
	m_year_birth $controls1 $controls2 i.FE_prov i.q_year if q_year!=1967 & q_year != 1968, cluster(prov)  nrtolerance(1e-2)
		
	* Build the Life Table
	keep if e(sample)
	keep if q_LLF==0
	
		* Replace LLF = 0 
		replace q_LLF = 0

		* Replace age*LLF interactions with zero and predict reference group
		forval i = 1/56 { 			
			replace LLF_quarter_`i' = 0
			} 
			
		cap drop xb
		predict xb, xb
		
		cap drop qx
		predict qx
		
			forval q = 1/56 { 
			
				sum xb if quarter == `q' 
				local xb = r(mean)
				
				sum qx if quarter == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_NoLLF = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_NoLLF = LT_NoLLF \ temp
					}
				}

		* Replace LLF = 1 
		replace q_LLF = 1

		* Replace age*LLF interactions with zero and predict LLF
		forval i = 1/56 { 			
			replace LLF_quarter_`i' = 1 if quarter_`i' == 1
			} 
			
		cap drop xb
		predict xb, xb
		
		cap drop qx
		predict qx
		
			forval q = 1/56 { 
			
				sum xb if quarter == `q' 
				local xb = r(mean)
				
				sum qx if quarter == `q' 
				local qx = r(mean)
				
				if `q' == 1 {
					mat LT_LLF = `q', `xb', `qx' 
					} 
				else if `q' > 1 { 
					mat temp = `q', `xb', `qx' 
					mat LT_LLF = LT_LLF \ temp
					}
				}

	* Compute Lx and Dx for each age with and without LLF
	mat temp = J(56, 1, .) 
	mat LT_NoLLF = LT_NoLLF, temp, temp
	mat LT_LLF = LT_LLF, temp, temp
		
	mat LT_NoLLF[1,4] = 1
	mat LT_LLF[1,4] = 1
	
	forval i = 1/56 {
		mat LT_NoLLF[`i',5] = LT_NoLLF[`i',3]*LT_NoLLF[`i',4]
		mat LT_LLF[`i',5] = LT_LLF[`i',3]*LT_LLF[`i',4]
		
		if `i' < 56 { 
			local j = `i' + 1
			mat LT_NoLLF[`j' , 4] = LT_NoLLF[`i' , 4] - LT_NoLLF[`i',5]
			mat LT_LLF[`j' , 4] = LT_LLF[`i' , 4] - LT_LLF[`i',5]
		} 
	}
		
	mat colnames LT_NoLLF = "Quarter" "XB_NoLLF" "Qx_NoLLF" "Lx_NoLLF" "Dx_NoLLF" 
	mat colnames LT_LLF = "Quarter" "XB_LLF" "Qx_LLF" "Lx_LLF" "Dx_LLF"
	
	mat temp = LT_LLF[1..., 2...]	
	mat colnames temp = "XB_LLF" "Qx_LLF" "Lx_LLF" "Dx_LLF"
	mat LifeTable_1 = LT_NoLLF, temp
	 
	
	* Plot the CDF
	clear 
	svmat LifeTable_1, names(col)
	gen Advanced_NoLLF = 1 - Lx_NoLLF
	gen Advanced_LLF = 1 - Lx_LLF
	
	twoway (connected Advanced_NoLLF Quarter, lpattern(solid) msize(zero)) ///
			(connected Advanced_LLF Quarter, lpattern(dash) msize(zero)) ///
			, legend(pos(6) order(1 "Before LLF" 2 "After LLF") rows(1)) ///
			 title("Marriage to First Birth") xtitle("Quarter at Risk") ytitle("Share of Population") ///
			 xlab(0(4)56) xline(14, lcolor(gs9)) ///
			 text(.2 14 "Target Birth Interval" "Under LLF", size(small) box fcolor(white) bcolor(white)) ylab(0(.1)1, format(%9.2f))
	
	graph export "$d_fig\Fig_A15a_BirthIntervals.jpg", replace  height(600) width(900)
	
	* Compute the change in mean age and median interval length (quarters) 
	local mean_NoLLF = 0
	local dx_cum_NoLLF = 0 
	local flag_NoLLF = 1

	forval i = 1/26 { 
		local j = `i' - 1
		local int = LT_NoLLF[`i', 1]*LT_NoLLF[`i', 5]
		local dx = LT_NoLLF[`i', 5]
		local mean_NoLLF = `mean_NoLLF' + `int' 
		local dx_cum_NoLLF = `dx_cum_NoLLF' + `dx' 
		
		local y = LT_NoLLF[`i', 4]
		if `y' < .5 & `flag_NoLLF' == 1 { 
			local flag_NoLLF = 2 
			local y_1 = LT_NoLLF[`j', 4]
			local x = LT_NoLLF[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_NoLLF = (.5-`alpha')/`beta'
			}
	}		
	
	local mean_LLF = 0 
	local dx_cum_LLF = 0 
	local flag_LLF = 1
	forval i = 1/28 { 
		local j = `i' - 1
		local int = LT_LLF[`i', 1]*LT_LLF[`i', 5]
		local dx = LT_LLF[`i', 5]
		local mean_LLF = `mean_LLF' + `int' 
		local dx_cum_LLF = `dx_cum_LLF' + `dx' 
		
		local y = LT_LLF[`i', 4]
		if `y' < .5 & `flag_LLF' == 1 { 
			local flag_LLF = 2 
			local y_1 = LT_LLF[`j', 4]
			local x = LT_LLF[`i', 1]
			local beta = `y' - `y_1'
			local alpha = `y'-(`beta'*`x')
			local median_LLF = (.5-`alpha')/`beta'
			}
		} 
		
	local mean_NoLLF = `mean_NoLLF'*`dx_cum_NoLLF'
	local mean_LLF = `mean_LLF'*`dx_cum_LLF'
	
	local dif_mean = `mean_LLF' - `mean_NoLLF'
	local dif_median = `median_LLF' - `median_NoLLF'
	
	mat effect1 = `mean_LLF', `mean_NoLLF', `dif_mean', `median_LLF', `median_NoLLF', `dif_median'
	mat colnames effect1 = "MeanLLF" "MeanNoLLF" "DiffMean" "MedianLLF" "MedianNoLLF" "DiffMedian"
	
	restore

	*********************************
	* (2.2) Second - Fourth births
	*********************************
	
	forval interval = 2/4 {
	
	preserve
		keep if order == `interval'  & year_marriage >= 1960 & birth_year <= 1979
		
		local prev = `interval' - 1

		keep marriage_age m_year_birth educ U5MR_5yr_ave $controls3 birth_month birth_year ///
		LLFyear year uniqueid interval`interval' interval`prev' prov marcmc FE_prov
		
		* I want to make this a woman-quarter level dataset so I can allow LLF to vary over the interval
		* gen quarter indicators that capture survival time - they should be zeros until the quarter the interval ends, when it should be 1
		gen 	quarter1=0 if interval`interval' >= 12 & interval`interval' != .
		replace quarter1=1 if interval`interval' < 12
		
		forval i = 2/28 { 
			local j 	= 3*`i' 
			local low 	= 9 + `j' - 3
			local high 	= 9 + `j'
		
			gen 	quarter`i'=0 if interval`interval'>=`high' & interval`interval' != .
			replace quarter`i'=1 if interval`interval'>=`low'  & interval`interval' < `high'	
		}
		
		
		* Now can identify the year and month at the start of the interval 
		* Note that first month at risk is 9 months after the previous event, but the start is the start of the entire birth interval
		gen end_cmc = 12*(birth_year-1900) + birth_month
		gen start_cmc = end_cmc - interval`interval' 
		
		forval i = 1/28 { 
			cap drop q_cmc`i'
			gen q_cmc`i' = .
			
			local q = `i'*3 - 3
			replace q_cmc`i' = start_cmc + 9 + `q' if quarter`i'  != .
			
			cap drop q_year`i'  
			gen q_year`i' = . 
			replace q_year`i' = q_cmc`i'/12 + 1900 if quarter`i'  != .
			
			cap drop q_LLF`i' 
			gen q_LLF`i' = . 
			replace q_LLF`i' = 0 if quarter`i' != . 
			replace q_LLF`i' = 1 if q_year`i' >= LLFyear & quarter`i' != . 
			} 
			
		drop q_cmc*  end_cmc start_cmc
		
		reshape long q_LLF q_year quarter, i(uniqueid) j(interval_q)
		replace q_year = floor(q_year)
		rename quarter birth
		rename interval_q quarter
		drop if birth == . 
		stset quarter, failure(birth) id(uniqueid)
			
		tab quarter, gen(quarter_)
		forval i = 1/28 { 
			gen LLF_quarter_`i' = q_LLF*quarter_`i' 
			} 
			
		gen birth_year2LLF = birth_year-LLFyear
		
		keep if q_year >= 1964
		keep if birth_year2LLF >= -4
		keep if birth_year2LLF <= 9		
		
		xi: logit 	birth q_LLF LLF_quarter_2-LLF_quarter_28 quarter_2-quarter_28 m_year_birth $controls1 $controls2 ///
					i.FE_prov i.q_year if q_year!=1967 & q_year != 1968, cluster(prov) iter(30) nrtolerance(1e-2)
			
		* Build the Life Table
		keep if e(sample)
		keep if q_LLF==0
		
			* Replace LLF = 0 
			replace q_LLF = 0

			* Replace age*LLF interactions with zero and predict reference group
			forval i = 1/28 { 			
				replace LLF_quarter_`i' = 0
				} 
				
			cap drop xb
			predict xb, xb
			
			cap drop qx
			predict qx
			
				forval q = 1/28 { 
				
					sum xb if quarter == `q' 
					local xb = r(mean)
					
					sum qx if quarter == `q' 
					local qx = r(mean)
					
					if `q' == 1 {
						mat LT_NoLLF = `q', `xb', `qx' 
						} 
					else if `q' > 1 { 
						mat temp = `q', `xb', `qx' 
						mat LT_NoLLF = LT_NoLLF \ temp
						}
					}

			* Replace LLF = 1 
			replace q_LLF = 1

			* Replace age*LLF interactions with zero and predict LLF
			forval i = 1/28 { 			
				replace LLF_quarter_`i' = 1 if quarter_`i' == 1
				} 
				
			cap drop xb
			predict xb, xb
			
			cap drop qx
			predict qx
			
				forval q = 1/28 { 
				
					sum xb if quarter == `q' 
					local xb = r(mean)
					
					sum qx if quarter == `q' 
					local qx = r(mean)
					
					if `q' == 1 {
						mat LT_LLF = `q', `xb', `qx' 
						} 
					else if `q' > 1 { 
						mat temp = `q', `xb', `qx' 
						mat LT_LLF = LT_LLF \ temp
						}
					}


		* Compute Lx and Dx for each age with and without LLF
		mat temp = J(28, 1, .) 
		mat LT_NoLLF = LT_NoLLF, temp, temp
		mat LT_LLF = LT_LLF, temp, temp
			
		mat LT_NoLLF[1,4] = 1
		mat LT_LLF[1,4] = 1
		
		forval i = 1/28 {
			mat LT_NoLLF[`i',5] = LT_NoLLF[`i',3]*LT_NoLLF[`i',4]
			mat LT_LLF[`i',5] = LT_LLF[`i',3]*LT_LLF[`i',4]
			
			if `i' < 28 { 
				local j = `i' + 1
				mat LT_NoLLF[`j' , 4] = LT_NoLLF[`i' , 4] - LT_NoLLF[`i',5]
				mat LT_LLF[`j' , 4] = LT_LLF[`i' , 4] - LT_LLF[`i',5]
			} 
		}
			
		mat colnames LT_NoLLF = "Quarter" "XB_NoLLF" "Qx_NoLLF" "Lx_NoLLF" "Dx_NoLLF" 
		mat colnames LT_LLF = "Quarter" "XB_LLF" "Qx_LLF" "Lx_LLF" "Dx_LLF"

		mat temp = LT_LLF[1..., 2...]	
		mat colnames temp = "XB_LLF" "Qx_LLF" "Lx_LLF" "Dx_LLF"	

		mat LifeTable_`interval' = LT_NoLLF, temp
		 
		
		* Plot the CDF	
		clear 
		svmat LifeTable_`interval', names(col)
		gen Advanced_NoLLF = 1 - Lx_NoLLF
		gen Advanced_LLF = 1 - Lx_LLF
		
		if `interval' == 2 { 
		twoway (connected Advanced_NoLLF Quarter, lpattern(solid) msize(zero)) ///
				(connected Advanced_LLF Quarter, lpattern(dash) msize(zero)) ///
				, legend(pos(6) order(1 "Before LLF" 2 "After LLF") rows(1)) ///
				 title("First to Second Birth Interval") xtitle("Quarter at Risk") ytitle("Share of Population") xlab(0(2)28) xline(14, lcolor(gs9)) ///
				 text(.2 14 "Target Birth Interval" "Under LLF", size(small) box fcolor(white) bcolor(white)) ylab(0(.1)1, format(%9.2f))
		graph export "$d_fig\Fig_A15b_BirthIntervals.jpg", replace  

		 }
		
		else if `interval' == 3 { 
		twoway (connected Advanced_NoLLF Quarter, lpattern(solid) msize(zero)) ///
				(connected Advanced_LLF Quarter, lpattern(dash) msize(zero)) ///
				, legend(pos(6) order(1 "Before LLF" 2 "After LLF") rows(1)) ///
				 title("Second to Third Birth Interval") xtitle("Quarter at Risk") ytitle("Share of Population") xlab(0(2)28) xline(14, lcolor(gs9)) ///
				 text(.2 14 "Target Birth Interval" "Under LLF", size(small) box fcolor(white) bcolor(white)) ylab(0(.1)1, format(%9.2f))
		graph export "$d_fig\Fig_A15c_BirthIntervals.jpg", replace  

		} 
		
		else if `interval' == 4 { 
		twoway (connected Advanced_NoLLF Quarter, lpattern(solid) msize(zero)) ///
				(connected Advanced_LLF Quarter, lpattern(dash) msize(zero)) ///
				, legend(pos(6) order(1 "Before LLF" 2 "After LLF") rows(1)) ///
				 title("Third to Fourth Birth Interval") xtitle("Quarter at Risk") ytitle("Share of Population") xlab(0(2)28) xline(14, lcolor(gs9)) ///
				 text(.2 14 "Target Birth Interval" "Under LLF", size(small) box fcolor(white) bcolor(white)) ylab(0(.1)1, format(%9.2f))
		graph export "$d_fig\Fig_A15d_BirthIntervals.jpg", replace  
		} 
			
		* Compute the change in mean age and median interval length (quarters) 
		local mean_NoLLF = 0
		local dx_cum_NoLLF = 0 
		local flag_NoLLF = 1

		forval i = 1/28 { 
			local j = `i' - 1
			local int = LT_NoLLF[`i', 1]*LT_NoLLF[`i', 5]
			local dx = LT_NoLLF[`i', 5]
			local mean_NoLLF = `mean_NoLLF' + `int' 
			local dx_cum_NoLLF = `dx_cum_NoLLF' + `dx' 
			
			local y = LT_NoLLF[`i', 4]
			if `y' < .5 & `flag_NoLLF' == 1 { 
				local flag_NoLLF = 2 
				local y_1 = LT_NoLLF[`j', 4]
				local x = LT_NoLLF[`i', 1]
				local beta = `y' - `y_1'
				local alpha = `y'-(`beta'*`x')
				local median_NoLLF = (.5-`alpha')/`beta'
				}
		}		
		
		local mean_LLF = 0 
		local dx_cum_LLF = 0 
		local flag_LLF = 1
		forval i = 1/28 { 
			local j = `i' - 1
			local int = LT_LLF[`i', 1]*LT_LLF[`i', 5]
			local dx = LT_LLF[`i', 5]
			local mean_LLF = `mean_LLF' + `int' 
			local dx_cum_LLF = `dx_cum_LLF' + `dx' 
			
			local y = LT_LLF[`i', 4]
			if `y' < .5 & `flag_LLF' == 1 { 
				local flag_LLF = 2 
				local y_1 = LT_LLF[`j', 4]
				local x = LT_LLF[`i', 1]
				local beta = `y' - `y_1'
				local alpha = `y'-(`beta'*`x')
				local median_LLF = (.5-`alpha')/`beta'
				}
			} 
			
		local mean_NoLLF = `mean_NoLLF'*`dx_cum_NoLLF'
		local mean_LLF = `mean_LLF'*`dx_cum_LLF'
		
		local dif_mean = `mean_LLF' - `mean_NoLLF'
		local dif_median = `median_LLF' - `median_NoLLF'
		
		mat effect`interval' = `mean_LLF', `mean_NoLLF', `dif_mean', `median_LLF', `median_NoLLF', `dif_median'
		mat colnames effect`interval' = "MeanLLF" "MeanNoLLF" "DiffMean" "MedianLLF" "MedianNoLLF" "DiffMedian"

	restore
}	
